The catenary is the curve formed by a homogeneous chain suspended between
two points in a uniform gravity field. It looks like a parabola, but it's
mathematically quite different. The equation of a parabola is just
 y = a*x*x, a simple quadratic, whereas a catenary is
 y = a*cosh(x/a), a transcendental equation.

I had occasion to look at Chris Colefax's chain building code in "Linc.inc",
after referring a new user to it, and was dismayed to discover that Chris
used quadratics to build his chains. :(

The Chain() macro below constructs a true catenary. This is a preliminary
draft, so there are no docs as yet, sorry. Read the comments for hints. :)

Any questions and comments are most welcome. Have fun!

//-------------------------------------------------------------------------

// Persistence of Vision Ray Tracer Include File
// File: Catenary.inc
// Vers: 3.6
// Desc: Proper catenary chain.
// Date: 2005.07.30
// Auth: PM 2Ring
//
// Catenary parameter calculations thanks to Zdislav V. Kovarik. See below
//
//-------------------------------------------------------------------------

#ifndef(Catenary_Inc_Temp)
#declare Catenary_Inc_Temp=version;
#version 3.5;

#ifdef(View_POV_Include_Stack)
    #debug "including Catenary.incn"
#end

//-------------------------------------------------------------------------
//
// From: kov### [at] mcmailcisMcMasterCA (Zdislav V. Kovarik)
// Subject: Re: Catenary
// Date: 5 Nov 1999 14:33:22 -0500
// Newsgroups: sci.math
// Keywords: fitting a catenary to match a suspended cable
//
// Upright catenary:  y - y_0 = a * cosh((x - x_0)/a) ,  a > 0.
// The vertex is (x_0, y_0+a).
// The parameter a turns out to be the radius of curvature at the vertex.
// Remark: The radius of curvature at a point above x is
//   a * (cosh((x-x_0)/a)^2.
//
// Problem: Parameters a, x_0, y_0 to be found so that the catenary arc
passes
// through (x_1, y_1), (x_2, y_2) and has length L between these points
// (provided L > sqrt((x_2 - x_1)^2 + (y_2 - y_1)^2) )
//
// The equations to be solved are (after some symmetrizing manipulation
// of the arclength integral)
//
//   y_1 - y_0 = a * cosh ((x_1 - x_0)/a)
//   y_2 - y_0 = a * cosh ((x_2 - x_0)/a)
//   2*a * cosh((x_1 + x_2 - 2*x_0)/(2*a)) * sinh((x_2 - x_1)/(2*a)) = L
//
// and it can be reduced to solving for one unknown at a time.
//
// First, we introduce an auxiliary unknown z which is to satisfy
//
//   sinh(z) / z = sqrt(L^2 - (y_2 - y_1)^2) / abs(x_2 - x_1)  ,   z > 0
//
// (the only transcendental non-elementary equation)
// and then the unknowns pop out:
//
//     a = abs(x_2 - x_1) / (2*z)
//   x_0 = (1/2)*(x_1 + x_2 - a * ln ((L + (y_2 - y_1)) / (L - (y_2 - y_1)))
//   y_0 =  (y_1 + y_2)/2 - (L/2) * coth(z).
//
//-------------------------------------------------------------------------

//Find B such that A=sinh(B)/B. Named in parallel to sinc()
#macro asinch(A)
  #local B = sqrt(6*max(1e-4,A-1));   //1st approx, from sinh(B) = B +
B^3/3! + ...
  #local B = asinh(A*B);              //2nd approx, from sinh(B) = A*B

  #local I = 0;
  #while(I<8)                         //Newton's method
    #local S = sinh(B);
    #local C = B * cosh(B);
    #local M = (A*B + C - 2*S) / (C - S);
    #local B = B * M;
    #local I = I + 1;
    #if(abs(M-1)<1e-12)               //bailout
      #local I=8;
    #end
  #end
  B                                   //return value
#end

//Find width of an object's bounding box
#macro BBWidth(A) (max_extent(A) - min_extent(A)).x #end

//Make a catenary. Parameters: Link object, Start point,End Point,
//Slackness of the chain >1, Link overlap, Extra twist on whole cable
#macro Chain(Link, StartA, EndA, Slack, Overlap, Twist)
  //Temporarily translate to origin & work in XY plane
  #local End = EndA - StartA;
  #local TH = degrees(atan2(End.z, End.x));
  #local End = vrotate(End, y*TH);              //Rotate to positive x-axis

  //Find required chain length and number of links
  #local Len = vlength(End) * Slack;            //Basic chain length
  #local LL = BBWidth(Link) / Overlap;          //Link Length adjusted for
link overlap
  #local Steps = 2*floor(.5*(.5+Len/LL));       //Round up to an even
integer number of links
  #local Len = Steps * LL;                      //Adjusted chain length

  //Find vertex of catenary that connects <0,0,0> & End, with length Len.
  #local P = sqrt(Len*Len - End.y*End.y) / End.x;
  #local Q = asinch(P);
  #local A = End.x / (2*Q);                     //Catenary curvature
parameter.
  #local X = A * ln((Len + End.y) / (Len - End.y));
  #local Y = Len/tanh(Q);
  #local V = (End - <X, Y, 0>)*.5;              //Vertex
  #local S1 = A*sinh(V.x/A);                    //Arclength at vertex

  //Step evenly along catenary parametrized by arclength.
  union {
    #local I=1;
    #while (I<Steps)
      #local S = LL * I - S1;                   //Arclength from Vertex
      #local M = S / A;                         //Slope of tangent
      #local X = asinh(M);
      #local Y = sqrt(1 + M*M);

      object{
        Link
        rotate (90*mod(I+1,2) + Twist*360*(I-1)/(Steps-2))*x
        rotate z*degrees(atan(M))               //Rotate parallel to tangent
        translate V + A*<X, Y, 0>
      }
      #local I=I+1;
    #end

    //Transform back
    rotate -y*TH
    translate StartA
  }
#end

#version Catenary_Inc_Temp;
#end

//-------------------------------------------------------------------------



//-------------------------------------------------------------------------

// Persistence of Vision Ray Tracer Scene Description File
// File: Catenary.pov
// Vers: 3.6
// Desc: Test Catenary include file
// Date: 2005.07.30
// Auth: PM 2Ring
//
// Catenary parameter calculations thanks to Zdislav V. Kovarik
// See Catenary.inc for details.
//
//-------------------------------------------------------------------------
//
// -f -A0.4 +AM2 +R1
// -d +A0.05 +AM2 +R3
//

#include "finish.inc"
#include "metals.inc"

//Chain making macro
#include "Catenary.inc"

global_settings {
  assumed_gamma 1.0
  max_trace_level 25
}

//-------------------------------------------------------------------------

//Simple chain macro. Parameters: Start point,End Point. Make sure other
items are declared before calling!
#macro ChainQ(Start, End) Chain(Link, Start, End, Slack, Overlap, Twist)
#end

//Chain terminal post
#macro Terminal(Pos)
union{
  sphere{Pos, PostRad*1.6}
  cylinder{Pos*<1,0,1>, Pos-0.35*PostRad*y, PostRad}

  pigment{rgb <.2, .5, 1>}
  finish{Glossy}
}
#end

//Chain, with terminal at start
#macro TermChain(Start, End)
  Terminal(Start)
  ChainQ(Start, End)
#end

//Link objects
#declare Torus = torus {.75, .175 scale 0.075*<1, 1, .65> }
#declare Torus1 = object {Torus scale 2 texture{T_Gold_2E} rotate 0*45*x}

//--- The scene -----------------------------------------------------------

#declare Rad = 2.0;               //Scene size control
#declare PostRad= 0.150;          //Post radius

//Chain parameters
#declare Link = Torus1;           //Link object
#declare Slack = 1.12;            //Slackness of the chain. (Length of
chain) / (straight distance between points)
#declare Overlap = 1.65;          //Link overlap
#declare Twist = 0;               //Chain twist (in cycles)

//Points to connect
#declare V1 = < 1.5*Rad, 0.75*Rad, 1>;
#declare V2 = <-1.5*Rad, 1.25*Rad, 3>;

//Do it!
TermChain(V1, V2)
Terminal(V2)

//Simple room with checkered floor
#declare WS = 5*Rad;
box{<-1, -2/WS, -2>, <1, 3, 2> scale WS inverse pigment{gradient y scale
3.001*WS}finish{Shiny}}
box{<-1, -1/WS, -2>, <1, 0, 2> scale WS pigment{checker rgb 1,rgb
....05}finish{Glossy diffuse 0.80}}

camera {
  location <-0.5, 3.5, -5.5> * 0.93 * Rad
  look_at  y*2

  right x*image_width/image_height up y
  direction z

  angle 30
}

light_source {<1, 9, -3>*Rad rgb 1 spotlight point_at z*2 falloff 16 radius
5 }

//-------------------------------------------------------------------------

Post a reply to this message

Attachments:
Download 'catenaryf2.jpg' (100 KB)

Preview of image 'catenaryf2.jpg'